\documentclass[a4paper]{article}
\usepackage[affil-it]{authblk}
% \usepackage[backend=bibtex,style=numeric]{biblatex}
\usepackage[UTF8]{ctex}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{geometry}
\usepackage{graphicx}
\usepackage{subcaption}
\usepackage{booktabs}
\usepackage{float}
\geometry{margin=1.5cm, vmargin={0pt,1cm}}
\setlength{\topmargin}{-1cm}
\setlength{\paperheight}{29.7cm}
\setlength{\textheight}{25.3cm}

\begin{document}
% ==================================================
\title{NA项目作业：报告}

\author{金郁香 3220103054
  \thanks{电子邮箱: \texttt{2621201771@qq.com}}}
\affil{浙江大学 信计2201班}

\date{截止时间: \today}

\maketitle

% ============================================
\section*{1. 在 ppForm 和 BSpline 格式下分别实现线性样条函数}

在Problem1.cpp中输入了knots = $\{0, 1, 2, 3, 4\}$和coefficients = $\{0, 2, 4, 6, 8\}$以测试两种线性样条对函数$y=2x$的插值效果。程序会输出样条函数在$x=0.5,1.5,2.5,3.5$处的插值结果，观察它们与理论值$y =1, 3, 5, 7$的差别，可知插值符合预期。

\begin{table}[H]
    \centering
    \caption{LinearSpline}
    \begin{tabular}{lll}
        \toprule 
        x & B-Spline & pp-Spline\\
        \midrule
        0.5 & 1 & 1\\
        1.5 & 3 & 3\\
        2.5 & 5 & 5\\
        3.5 & 7 & 7\\
        \bottomrule
    \end{tabular}
\end{table}

\section*{2. 在 ppForm 格式下分别推导并实现三类三次样条函数}

在Problem2.cpp中输入了knots = $\{0, 0.3, 1, 1.6, 2\}$和coefficients = $sin(\pi*\{knots\})$以测试pp形式下三类三次样条函数对函数$y=sin(\pi x)$的插值效果。程序会输出样条函数在$x=\frac{1}{6}, \frac{1}{4}, \frac{1}{3}$处的插值结果，观察它们与理论值$y =0.5, 0.707, 0.866$的差别，可知插值基本符合预期。

\begin{table}[H]
    \centering
    \caption{pp-3-Spline}
    \begin{tabular}{llll}
        \toprule 
        x & Periodic & Natural & Clamped\\
        \midrule
        0.166667 & 0.512971 & 0.514539 & 0.508772\\
        0.25 & 0.716698 & 0.717329 & 0.714962\\
        0.333333 & 0.855415 & 0.855069 & 0.856426\\
        \bottomrule
    \end{tabular}
\end{table}

\section*{3. 在BSpline 格式下分别推导并实现如下三类三次样条函数 }

在Problem3.cpp中输入了knots = $\{0, 0.3, 1, 1.6, 2\}$和coefficients = $sin(\pi*\{knots\})$以测试三类三次B样条函数对函数$y=sin(\pi x)$的插值效果。程序会输出样条函数在$x=\frac{1}{6}, \frac{1}{4}, \frac{1}{3}$处的插值结果，观察它们与理论值$y =0.5, 0.707, 0.866$的差别，可知插值基本符合预期。

\begin{table}[H]
    \centering
    \caption{B-3-Spline}
    \begin{tabular}{llll}
        \toprule 
        x & Periodic & Natural & Clamped\\
        \midrule
        0.166667 & 0.512971 & 0.514539 & 0.508772\\
        0.25 & 0.716698 & 0.717329 & 0.714962\\
        0.333333 & 0.855415 & 0.855069 & 0.856426\\
        \bottomrule
    \end{tabular}
\end{table}

\section*{4. 验证 ppForm 和 BSpline 格式下取相同插值点和相同边界条件，得到的曲线是相同的}

观察在Problem2、3（相同插值点和相同边界条件）中两种样条函数在相同点处的插值结果，已经可以注意到它们是完全相同的。为了进一步验证结果，在Problem4.cpp中取了一组更加随意的数据分别输入给 ppForm 和 BSpline 格式下的三类样条函数，并输出两者在相同点处的插值结果。观察输出结果发现两者完全相同，这表明两种样条插值得到的曲线是一样的。

\begin{table}[H]
    \centering
    \caption{pp-3-Spline Vs B-3-Spline}
    \begin{tabular}{lllllll}
        \toprule 
        x & pp-Periodic & B-Periodic& pp-Natural & B-Natural & pp-Clamped & B-Clamped\\
        \midrule
        0.2 & 1.23409 & 1.23409 & 1.23389 & 1.23389 & 1.21983 & 1.21983\\
        1.35 & 0.95328 & 0.95328 & 0.938892 & 0.938892 & 0.919454 & 0.919454\\
        2 & 0.203813 & 0.203813 & 0.284156 & 0.284156 & 0.409161 & 0.409161\\
        3.7 & 0.626505 & 0.626505 & -0.58825 & -0.58825 & -1.24919 & -1.24919\\
        \bottomrule
    \end{tabular}
\end{table}

\section*{5. 实现 BSpline 格式在任意阶、任意节点样条的绘制。}

在Problem5.cpp中，输入knots = $\{1, 2, 3, 4, 5, 6, 7\}$和coefficients = $\{1, 1, 1, 1, 1, 1, 1, 1, 1, 1\}$，以测试任意阶数任意节点的B样条能否绘制出理想的曲线：$n = 4, \, S(t) = \sum_{i=-2}^{7}1\cdot B_i^4(t) == 1, \, t\in[1,7]$。程序会输出函数在$x=1.3, 1.5, 2.4, 3.5, 4.1, 5.8, 6.9$处的函数值，观察它们与理论值$y =1,1,1,1,1,1,1$的差别，可知该B样条求和公式基本符合预期。

\begin{table}[H]
    \centering
    \caption{LinearSpline}
    \begin{tabular}{lc}
        \toprule 
        x & Sum of B-Spline functions\\
        \midrule
        1.3 & 1 \\
        1.5 & 1 \\
        2.4 & 1 \\
        3.5 & 1 \\
        4.1 & 1 \\
        5.8 & 1 \\
        6.9 & 1 \\
        \bottomrule
    \end{tabular}
\end{table}

另外，还利用该公式绘制出了如下0到4阶B样条基函数图像：

\begin{figure}[H]
    \centering
    \includegraphics[width=0.4\textwidth]{figure/P5_1.png}
    \caption{x=3处的0到4阶B样条基函数图像}
\end{figure}

下面再绘制一个给定不均匀节点和对应系数的$S(t)$图像，这里取n = 5, knots = $\{1, 1.6, 3, 4.2, 5.3, 6.6, 7\}$和coefficients = $\{3.2, 1, 5, 2, -1, 0.1, 3, 4, 1, -8, 3\}$，得到结果如下图，可知该B样条求和公式基本符合预期。

\begin{figure}[H]
    \centering
    \includegraphics[width=0.4\textwidth]{figure/P5_2.png}
    \caption{不均匀节点下的B样条求和函数图像}
\end{figure}

\section*{6. 实现平面上的样条曲线拟合}

测试结果详情见E题的r1、r2部分。

\section*{7. 实现球面上的样条曲线拟合}

测试结果详情见E题的r3部分。

在此给出球面上样条曲线拟合的思路说明：球面上的曲线虽然由x,y,z三个坐标决定，但x,y,z又可以通过球面参数方程由u,v两个变量来表示。所以我们只需要对u,v进行平面曲线拟合，再用得到的u(t),v(t)来获取x(t),y(t),z(t)的值，从而实现球面曲线的拟合。

其中，用于坐标变换的球面参数方程（R给定）如下；
\[ \left\{ \begin{array}{l}
x = R\sin u \cos v, \\
y = R\sin u \sin v, \\
z = R\cos u,
\end{array} \right. \]

\section*{A. 使用三次样条插值拟合函数}

在ProblemA.cpp中，这里选用了pp形式的完全三次样条在区间 $[-1,1]$ 上插值函数 $f(x) = \frac{1}{1+25x^2}$，并使用 $N=6,11,21,41,81$ 来观察不同的节点数下的插值效果。同时计算了每个子区间中点处的样条插值误差，并输出最大误差如下：

\begin{table}[H]
    \centering
    \caption{Max-norm Errors}
    \begin{tabular}{ll}
        \toprule 
        N & Max-norm error \\
        \midrule
        6 & 0.421705\\
        11 & 0.0205289 \\
        21 & 0.00316894 \\
        41 & 0.000275356 \\
        81 & 1.609e-05 \\
        \bottomrule
    \end{tabular}
\end{table}

结合上述结果和以下样条插值函数和真实函数的图像发现，除了$N = 6$时节点数量太少导致插值误差较大，在其余$N$的取值下，三次样条插值的拟合效果都较好，并未出现荣格现象；且随着节点数的增加，拟合准确度可以达到很高。

% 插入图片
\begin{figure}[H]
    \centering
    \includegraphics[width=0.4\textwidth]{figure/PA_6.png}
    \caption{N=6 时的样条插值}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.4\textwidth]{figure/PA_11.png}
    \caption{N=11 时的样条插值}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.4\textwidth]{figure/PA_21.png}
    \caption{N=21 时的样条插值}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.4\textwidth]{figure/PA_41.png}
    \caption{N=41 时的样条插值}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.4\textwidth]{figure/PA_81.png}
    \caption{N=81 时的样条插值}
\end{figure}
% \clearpage

\section*{C. 分别用二次和三次B样条进行函数插值}

这里分别用二次B样条函数的和三次B样条函数在区间 $[-5,5]$ 上对函数 $f(x) = \frac{1}{1+x^2}$进行插值。观察插值结果发现，在区间内的大部分区域（除了$x=0$附近），两个样条插值函数的拟合效果都较好，其中三次样条的插值效果要优于二次样条。

两个样条插值函数和真实函数的对比图如下：

% 插入图片
\begin{figure}[H]
    \centering
    \includegraphics[width=0.4\textwidth]{figure/PC.png}
    \caption{2、3次B样条插值函数与真实函数的对比图}
\end{figure}

\section*{D. 分析两种B样条插值函数的误差}

利用程序输出误差结果如下：

\begin{table}[H] % 设置表格浮动位置为“此处”（here）
\centering % 表格居中显示
\caption{Interpolation Errors} % 表格标题
\begin{tabular}{lll} % 定义表格列格式，l为左对齐，r为右对齐，r为居中对齐
    \toprule % 表格顶部横线
    x & Cubic & Quadratic \\ % 表头
    \midrule % 表格中部横线
    -3.5 & 0.000669568 & 0 \\
    -3   & 0 & 0.00141838 \\
    -0.5 & 0.0205289 & 1.11022e-16 \\
    0   & 0 & 0.120238 \\
    0.5 & 0.0205289 & 1.11022e-16 \\
    3   & 1.38778e-17 & 0.00141838 \\
    3.5 & 0.000669568 & 0 \\
    \bottomrule % 表格底部横线
\end{tabular}
\end{table}

观察发现，两种样条插值函数在它们各自的插值节点处——其中三次样条为整数点：$x = -3,0,3$，二次样条为：$x = -3.5,-0.5,0.5,3.5$，插值误差接近于或小于机器精度。这表明在这些点上，样条插值在数学意义上的确等于真实函数值。

比较两种样条插值在这些点上的误差，我们发现在总体上，三次样条插值的误差更小，且二次样条插值在$x=0$处误差特别大。综上，可以判断三次B样条的插值结果更为精确。

\section*{E. 实现样条曲线拟合并对不同方法进行比较}
\subsection*{E-a.心形曲线拟合}
首先实现课本上的要求，对心形曲线分别采用两种参数的样条曲线拟合和贝塞尔曲线拟合。其中，由于曲线是闭合的，故在进行三次样条插值时选用的是周期边界条件。

最终得到在$N = 10,40,160$时的三种曲线拟合结果如下：

% 插入图片
\begin{figure}[H]
    \centering
    \includegraphics[width=0.4\textwidth]{figure/PE_1_10.png}
    \caption{N=10时的心形曲线拟合}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.4\textwidth]{figure/PE_1_40.png}
    \caption{N=40时的心形曲线拟合}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.4\textwidth]{figure/PE_1_160.png}
    \caption{N=160时的心形曲线拟合}
\end{figure}

观察到样条曲线和贝塞尔曲线的区别：

1.贝塞尔曲线是先将整个曲线分成了N段，再在每段区间上分别拟合小曲线，最后进行连接，因此在N较小时，拟合效果不佳。在每段区间上，由两个端点的函数值和导数值，产生了四个控制点并生成三次函数来拟合曲线，由于控制点受导数影响较大，且导数在定义域上不一定存在（可能是无穷），所以贝塞尔曲线拟合存在一定局限性。

2.样条曲线是在整个曲线上取出节点，一次性实现对整段曲线的拟合，拟合出来的曲线在整体上更光滑，也更符合真实曲线。且样条曲线只需要利用节点处的函数值，不需要获取额外的导数信息，因此适用范围更广。

3.观察到当N较小时，两种曲线在心形的尖角处拟合效果都不是很好，不过这是由心形曲线的本身的几何性质决定的（尖角处的形状与三次多项式相差过大）。

\subsection*{E-b.两种参数下的样条曲线拟合}
下面对项目作业额外要求的另外两条曲线在$N = 10,40,160$时分别用两种参数的样条插值曲线进行拟合，得到结果如下：

\subsubsection*{1. r2曲线拟合}

由于r2曲线不再是封闭的，这里使用了完全边界条件下（手工计算了两个端点处的导数值并输入）的三次B样条插值来拟合曲线。

\begin{figure}[H]
    \centering
    \includegraphics[width=0.4\textwidth]{figure/PE_2_10.png}
    \caption{N=10时的r2曲线拟合}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.4\textwidth]{figure/PE_2_40.png}
    \caption{N=40时的r2曲线拟合}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.4\textwidth]{figure/PE_2_160.png}
    \caption{N=160时的r2曲线拟合}
\end{figure}

观察发现，用累积弦长s作为参数拟合出来的曲线在原点附近不太光滑，这也许是因为原点附近，在相同的弦长变化下角度变化会很大（相当于插值的相邻节点间的距离过大），从而导致拟合误差较大；而在其余累积弦长不那么小的区域，两种参数下的样条曲线拟合效果接近，其中s参数下的样条曲线相对更接近真实曲线。

除了$N = 10$时插值节点数过少导致误差偏大，其余情况下，两种方法都产生了较好的曲线拟合效果。

\subsubsection*{2. r3曲线拟合}

由于三维曲线r3仍然是封闭的，这里使用周期边界条件进行三次B样条插值来拟合曲线。

\begin{figure}[H]
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
        \includegraphics[width=\textwidth]{figure/PE_3_10_t.png}
        \caption{N=10时t参数下的r3曲线拟合}
    \end{subfigure}
    \begin{subfigure}[b]{0.4\textwidth}
        \includegraphics[width=\textwidth]{figure/PE_3_10_s.png}
        \caption{N=10时s参数下的r3曲线拟合}
    \end{subfigure}
    \caption{N=10}
\end{figure}

\begin{figure}[H]
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
        \includegraphics[width=\textwidth]{figure/PE_3_40_t.png}
        \caption{N=40时t参数下的r3曲线拟合}
    \end{subfigure}
    \begin{subfigure}[b]{0.4\textwidth}
        \includegraphics[width=\textwidth]{figure/PE_3_40_s.png}
        \caption{N=40时s参数下的r3曲线拟合}
    \end{subfigure}
    \caption{N=40}
\end{figure}

\begin{figure}[H]
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
        \includegraphics[width=\textwidth]{figure/PE_3_160_t.png}
        \caption{N=160时t参数下的r3曲线拟合}
    \end{subfigure}
    \begin{subfigure}[b]{0.4\textwidth}
        \includegraphics[width=\textwidth]{figure/PE_3_160_s.png}
        \caption{N=160时s参数下的r3曲线拟合}
    \end{subfigure}
    \caption{N=160}
\end{figure}

观察结果发现，两种参数下的样条曲线拟合在三个N取值下拟合效果都很好。

\section*{F. 利用截断函数的差商证明B样条基函数和截断函数的关系}

通过 ProblemF.cpp 分别绘制了以下$n = 1,2$, $x \in [0,4]$时的截断函数图像，这里为了方便对照，绘图时以均匀节点为例。直接观察最后一列的差商函数，乘上$(t_{i+n}-t_{i-1})$即$(n+1)$后，与前面题目绘制过的B样条基函数对比，可知定理3.32成立。

\subsection*{1. n = 1}
\begin{figure}[H]
    \centering
    \begin{subfigure}{0.3\textwidth}
        \includegraphics[width=\linewidth]{figure/PF_1_11.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.3\textwidth}
        \includegraphics[width=\linewidth]{figure/blank.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.3\textwidth}
        \includegraphics[width=\linewidth]{figure/blank.png}
    \end{subfigure}
    \par\bigskip
    \begin{subfigure}{0.3\textwidth}
        \includegraphics[width=\linewidth]{figure/PF_1_12.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.3\textwidth}
        \includegraphics[width=\linewidth]{figure/PF_1_21.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.3\textwidth}
        \includegraphics[width=\linewidth]{figure/blank.png}
    \end{subfigure}
    \par\bigskip
    \begin{subfigure}{0.3\textwidth}
        \includegraphics[width=\linewidth]{figure/PF_1_13.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.3\textwidth}
        \includegraphics[width=\linewidth]{figure/PF_1_22.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.3\textwidth}
        \includegraphics[width=\linewidth]{figure/PF_1_31.png}
    \end{subfigure}
    \caption{Divided diﬀerence of truncated power functions of n = 1}
\end{figure}

\subsection*{2. n = 2}
\begin{figure}[H]
    \centering
    \begin{subfigure}{0.2\textwidth} % 调整宽度以适应4列
        \includegraphics[width=\linewidth]{figure/PF_2_11.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.2\textwidth}
        \includegraphics[width=\linewidth]{figure/blank.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.2\textwidth}
        \includegraphics[width=\linewidth]{figure/blank.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.2\textwidth}
        \includegraphics[width=\linewidth]{figure/blank.png}
    \end{subfigure}
    
    \par\bigskip
    \begin{subfigure}{0.2\textwidth}
        \includegraphics[width=\linewidth]{figure/PF_2_12.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.2\textwidth}
        \includegraphics[width=\linewidth]{figure/PF_2_21.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.2\textwidth}
        \includegraphics[width=\linewidth]{figure/blank.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.2\textwidth}
        \includegraphics[width=\linewidth]{figure/blank.png}
    \end{subfigure}
    
    \par\bigskip
    \begin{subfigure}{0.2\textwidth}
        \includegraphics[width=\linewidth]{figure/PF_2_13.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.2\textwidth}
        \includegraphics[width=\linewidth]{figure/PF_2_22.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.2\textwidth}
        \includegraphics[width=\linewidth]{figure/PF_2_31.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.2\textwidth}
        \includegraphics[width=\linewidth]{figure/blank.png}
    \end{subfigure}
    
    \par\bigskip
    \begin{subfigure}{0.2\textwidth}
        \includegraphics[width=\linewidth]{figure/PF_2_14.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.2\textwidth}
        \includegraphics[width=\linewidth]{figure/PF_2_23.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.2\textwidth}
        \includegraphics[width=\linewidth]{figure/PF_2_32.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.2\textwidth}
        \includegraphics[width=\linewidth]{figure/PF_2_41.png}
    \end{subfigure}
    
    \caption{Divided difference of truncated power functions of n = 2}
\end{figure}
% ===============================================

\end{document}